LAB4

1.1

downsampled:

upsampled:

import sys
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import numpy as np
import timeit
import cv2 as cv 
from skimage import exposure
from scipy.fftpack import fft, fftshift, ifft, fft2, ifft2 # can replace scipy.fft to older version called scipy.fftpack
import scipy
from scipy import signal
from scipy.signal import cont2discrete, lti, dlti, dstep
from numpy import cos,pi,real
from matplotlib.pyplot import subplot, figure, title, stem, axes, plot, imshow
from scipy.signal import   butter , lfilter, freqz
from scipy.interpolate import CubicSpline

def MSE_SNR_PSNR(SignalA, SignalB):
    
    mse = (np.square(SignalA - SignalB)).mean()
    
    error_ = (SignalA - SignalB)
    sig_mean = (np.square(SignalA)).mean()
    noise_mean = (np.square(error_)).mean()

    snr = 10*np.log10(sig_mean/noise_mean)
    psnr = 20*np.log10(255/np.sqrt(mse))  # 255 is the maximum grey value

    return mse, snr, psnr


def convert(img, target_type_min, target_type_max, target_type):
    imin = img.min();
    imax = img.max();

    a = (target_type_max - target_type_min) / (imax - imin);
    b = target_type_max - a * imax;
    new_img = (a * img + b).astype(target_type);
    return new_img

def downsample(sig, factor, Anti_Aliasing=False, img=False):
    if Anti_Aliasing:
        if img:
            lpf = np.ones((3,3))/9
            sig = cv.filter2D(sig, -1, lpf)
        else:
            b, a = butter(15, 1/factor, btype='low', analog=False)
            sig = lfilter(b, a, sig) 
    if img:
        sig = sig[::factor, ::factor]
    else:
        sig = sig[::factor]
    return sig

def uint8c(x):
    """custom uint8 - (if val < 0 then val = 0)"""
    x2 = np.copy(x)
    x2  = x2.clip(0, 255)
    return np.uint8(x2)

def sample_test(x,factor, printout=True):
    print(f"")
    x_downsampled = downsample(x, factor)
    x_downsampled_antialiasing = downsample(x, factor, True)
    x_fft = np.abs(fft(x))
    x_downsampled_fft = np.abs(fft(x_downsampled))
    x_downsampled_antialiasing_fft = np.abs(fft(x_downsampled_antialiasing))
    f1 = np.arange(-len(x_fft)/2, len(x_fft)/2)
    f2 = np.arange(-len(x_downsampled_fft)/2, len(x_downsampled_fft)/2)
    f3 = np.arange(-len(x_downsampled_antialiasing_fft)/2, len(x_downsampled_antialiasing_fft)/2)

    if printout:
        plt.figure(figsize=(15,10))
        subplot(321); stem(x); title('original')
        subplot(322); stem(f1, fftshift(x_fft)); title('original fft')
        
        subplot(323); stem(x_downsampled); title(f'decimation factor = {factor}')
        subplot(324); stem(f2, fftshift(x_downsampled_fft)); title('fft decimation')
        
        subplot(325);stem(x_downsampled_antialiasing);title('decimation antialiasing')
        subplot(326); stem(f3, fftshift(x_downsampled_antialiasing_fft)); title('decimation antialiasing fft')
        plt.tight_layout()
    return x_downsampled, x_downsampled_antialiasing
    

Signal sin - Acos(2pi10t) + 2Acos(2pi17t)+3Acos(2pi5t)

time = 1 # sec
fs = 200 # Hz
f = 10 # Hz
A = 7
Ts = 1/fs
t = np.arange(0, time-Ts, Ts)    
x = A*cos(2*pi*10*t) + 2*A*cos(2*pi*17*t)+3*A*cos(2*pi*5*t);
sample_test(x, 2)
sample_test(x, 4)
sample_test(x, 8)
sample_test(x, 17)

(array([ 42.        , -27.76897101,  20.96755674, -17.0690992 ,
         -9.52913744,  17.03932492, -21.05361932,  37.90040574,
        -27.84310872,  19.41239557, -20.31479123,  -2.58612709]),
 array([ 5.52411930e-15,  2.57472125e-04,  3.39310392e-01,  4.89487863e+00,
        -3.31853335e+00, -6.24926302e+00,  1.44132060e+01, -1.89891894e+01,
         1.84453961e+01, -1.38034405e+01,  5.80662156e+00,  3.68393273e+00]))

Row from image down sample 2,4,8,17

img = cv.imread(r"C:\Users\rom21\OneDrive\Desktop\imageLAB4\cameraman.tif",0)
x=img[70,:]
sample_test(x, 2)
sample_test(x, 4)
sample_test(x, 8)
sample_test(x, 17)

(array([162, 168, 171, 177,  14,  16,  79, 144,  63, 134, 168, 169, 163,
        157, 151, 144], dtype=uint8),
 array([ 2.13073173e-14,  1.12114420e-03,  2.08289626e+00,  5.86659563e+01,
         1.90954452e+02,  1.66223053e+02,  1.74678454e+02,  9.83228354e+01,
        -1.99445170e+01,  6.29735645e+01,  1.10052834e+02,  1.13039900e+02,
         1.05583938e+02,  1.91917891e+02,  1.68573857e+02,  1.66273204e+02]))

My signal can see aliasing

x = A*cos(2*pi*10*t) + 2*A*cos(2*pi*20*t) + 2*A*cos(2*pi*30*t)  
sample_test(x, 2)
sample_test(x, 4)
sample_test(x, 8)
sample_test(x, 17)

(array([ 35.        , -13.52653238,  -2.16311896,  -3.56015122,
          2.98935688, -14.        ,   5.66311896,  26.21262707,
        -20.48935688,   4.87405654,  -7.        ,   4.87405654]),
 array([ 4.60343275e-15,  1.62427404e-04,  1.37472032e-01,  1.19119573e+00,
        -3.18374634e-01, -2.94432192e-02,  1.19656826e-01, -8.96916238e-02,
         7.07075329e-02, -6.54555921e-02,  5.92103586e-02, -3.74673268e-02]))

Upsample ZOH, FOH, linear, cubic and sinc

def upsample1d(sig, factor, method = 'ZeroOrderHold'):
    if method == 'ZeroOrderHold':
        sig_upsampled = np.repeat(sig, factor)

    elif (method == 'FirstOrderHold'):
        sig_upsampled = np.zeros(len(sig) * factor) 
        sig_upsampled[::factor] = sig
        filter_linear = np.zeros(factor*2 -1)
        for i in range(int((len(filter_linear) + 1 )/2)):
            filter_linear[i] = filter_linear[int(len(filter_linear)-1-i)]  = (i+1)/factor
             #= filter_linear[int(len(filter_linear)-1 -i)]
        sig_upsampled = np.convolve(sig_upsampled, filter_linear, 'same') 
    
    elif (method == 'sinc'):
        sig_upsampled = np.zeros(len(sig) * factor) 
        sig_upsampled[::factor] = sig
        l = 50
        x = np.linspace(-l, l, factor*l-1)
        sinc_filter =np.sinc((np.pi/(factor+2.1))*x)
        sig_upsampled = np.convolve(sig_upsampled, sinc_filter, 'same')
        
    elif (method == 'cubic'):
        sig_upsampled = np.zeros(len(sig) * factor) 
        for i in range(0,len(sig)-3):
            f = CubicSpline([0, 1, 2], sig[i:i+3], bc_type='natural')
            x_new = np.linspace(0, 2, factor+1)
            sig_upsampled[i*factor:i*factor+factor+1] = f(x_new)
    return sig_upsampled

x = A*cos(2*pi*10*t) + 2*A*cos(2*pi*17*t)+3*A*cos(2*pi*5*t);

factor = 4
x_downsampled, x_downsampled_antialiasing = sample_test(x, factor, False)
#s3 = upsample1d(x_downsampled, factor, method = 'ZeroOrderHold')
#s4 = upsample1d(x_downsampled_antialiasing, factor, method = 'FirstOrderHold')
#s5 = upsample1d(x_downsampled_antialiasing, factor, method = 'Ideal')
methods= ['ZeroOrderHold', 'FirstOrderHold', 'sinc', 'cubic']
for method in methods:
    upsampled = upsample1d(x_downsampled, factor, method = method)
    upsampled = upsampled[:len(upsampled)-1]
    MSE = np.sum(np.square(x - upsampled) / len(x))
    SNR = 10*np.log10(np.std(x) / np.std(x - upsampled))
    fig = plt.figure()
    axes = fig.subplots(nrows=2, ncols=1)
    axes[0].plot(range(0, len(x)), x, range(0, len(x)), upsampled)
    axes[0].set(title = 'original signal (blue), upsampled ' + method + ' (orange)\nMSE = ' + "{:.3f}".format(MSE) + ', SNR = ' + "{:.2f}".format(SNR) + 'dB')
    axes[1].plot(x_downsampled)
    axes[1].set(title = 'downsampled signal')
    fig.tight_layout()

1.2

Implament quantize:

def quantize(signal, k):
    temp = signal
    if type(signal[0]) == np.float64:
        temp = uint8c(signal*255) 
    delta = 256/k
    signal_q =  uint8c(np.int64(np.double(temp) / delta) * delta + delta/2)

    if type(signal[0]) == np.float64:
        signal_q = np.double(signal_q / 255)
    return signal_q


def quatize_analaze(sig, num_bits):
    sig_float_q = quantize(sig, num_bits)
    err = sig - sig_float_q
    avg_err = np.sqrt(np.average(np.square(err)))
    snr_db = 20*np.log10(np.std(sig) / np.std(err))
    fig, axi = plt.subplots(2,1)
    axi[0].hist(sig.flatten(), bins = 255)
    axi[0].hist(sig_float_q.flatten(), bins = 255)
    axi[0].set(title = 'Histogram')
    axi[0].legend(['Float64','Float64 Quantized ' + str(num_bits) + '-bit'])
    axi[1].plot(np.arange(0, len(sig)), sig,np.arange(0,len(sig_float_q)), sig_float_q)
    axi[1].set(title = 'Signals\nSNR = ' + "{:.1f}".format(snr_db) + 'dB, Avg Error = ' + "{:.4f}".format(avg_err))
    axi[1].legend(['Float64','Float64 Quantized ' + str(num_bits) + '-bit'])
    fig.tight_layout()
    plt.show()
    return sig_float_q

All plot: uint sin 256 levels, double sin 256 levels, row from image and (uint sin 256 levels) /4

fs = 200
L = 1
f = 2
Ts = 1/fs
t = np.arange(0,L-Ts,Ts)
y = np.cos(2*np.pi*f*t)
y = uint8c((1 + y) * 128)
Min = np.min(y)
Max = np.max(y)
y_double = ((np.double(y) - Min) / (Max - Min))
img=cv.imread(r"C:\Users\rom21\OneDrive\Desktop\imageLAB4\cameraman.tif", 0)
factor=17
x=img[70,:]
y_4 = y/4 +96
num_bits = 4
quatize_analaze(x, num_bits)
quatize_analaze(y_double, num_bits)
quatize_analaze(y, num_bits)
quatize_analaze(uint8c(y_4), num_bits)
target_noise_db = 10
target_noise_watts = 10 ** (target_noise_db / 10)
mean_noise = 0
noise = np.random.normal(mean_noise, np.sqrt(target_noise_watts), len(y))
noise_sig = y+noise
noise_sig =((np.double(noise_sig) - Min) / (Max - Min))
sig_float_q = quatize_analaze(noise_sig, 2)

2. Decimation and iterpolation:

2.1.3

a. Downsample factors = [2, 4, 10]:

img = cv.imread(r"C:\Users\rom21\OneDrive\Desktop\imageLAB4\chrono.tif", 0)
chrono = img[1601:2500, 1201:2300]
height, width = chrono.shape
factors = [2, 4, 10]
down_chrono_list = []
#downsampled
plt.figure(figsize=(15,15))
plt.title(f"orig")
plt.imshow(chrono, cmap = 'gray', vmin = 0, vmax = 255, aspect='equal') 

for i, factor in enumerate(factors):
        downsampled_img = downsample(chrono, factor, img=True)
        plt.figure(figsize=(15,15))
        plt.title(f"factor:{factor}")
        plt.imshow(downsampled_img, cmap = 'gray', vmin = 0, vmax = 255,aspect='equal') 
        
        

For a 1: 2 aspect ratio it will be difficult to see the differences but it will be enough to take a look at the size of the matrix that represents the image to understand that we have shrunk it while losing information along the way of course. For decimation in a ratio of 1: 4 and of course for decimation in a ratio of 1:10 we can see the effects of matrix contraction.

b-c. iterpolation bilinear and ZOH for image

def Zero_Order_Hold2d(img, factor):
    new_img = np.zeros(tuple(np.array(img.shape)*factor))
    for c in range(img.shape[0]):
        new_img[c*factor:c*factor+factor,:] = np.repeat(img[c,:], factor)
    return new_img

def bilinear(image, factor):
  img_height, img_width = image.shape
  height, width = tuple(np.array(image.shape)*factor)
  image = image.ravel()
  x_ratio = float(img_width - 1) / (width - 1) 
  y_ratio = float(img_height - 1) / (height - 1) 
  y, x = np.divmod(np.arange(height * width), width)
  x_l = np.floor(x_ratio * x).astype('int32')
  y_l = np.floor(y_ratio * y).astype('int32')
  x_h = np.ceil(x_ratio * x).astype('int32')
  y_h = np.ceil(y_ratio * y).astype('int32')
  x_weight = (x_ratio * x) - x_l
  y_weight = (y_ratio * y) - y_l
  a = image[y_l * img_width + x_l]
  b = image[y_l * img_width + x_h]
  c = image[y_h * img_width + x_l]
  d = image[y_h * img_width + x_h]
  resized = a * (1 - x_weight) * (1 - y_weight) + \
            b * x_weight * (1 - y_weight) + \
            c * y_weight * (1 - x_weight) + \
            d * x_weight * y_weight
  return resized.reshape(height, width)

#%% 2.1.3 Interpolation Zero_Order_Hold 2d

for i, factor in enumerate(factors):
    downsampled_img = downsample(chrono, factor, img=True)
    resized = Zero_Order_Hold2d(downsampled_img, factor)
    #resized = cv.resize(downsampled_img, img.shape, interpolation=cv.INTER_LINEAR)
    plt.figure(figsize=(15,15))
    plt.title(f"factor {factor}, Zero_Order_Hold2d Interpolation")
    plt.imshow(resized, cmap = 'gray', vmin = 0, vmax = 255,aspect='equal')

#%% 2.1.3 Interpolation Bilinear 2d
dpi = [1250, 300, 150, 72]
for i, factor in enumerate(factors):
    downsampled_img = downsample(chrono, factor, img=True)
    resized = bilinear(downsampled_img, factor)
    plt.figure(figsize=(15,15))
    plt.title(f" Bilinear Interpolation {factor}")
    plt.imshow(resized, cmap = 'gray', vmin = 0, vmax = 255,aspect='equal')

2.2.1 quantization quant_levels = [2, 4, 8, 16, 32, 64, 128] to skull image

        
image_path = r"C:\Users\rom21\OneDrive\Desktop\imageLAB4\skull.JPG"

img = cv.imread(image_path)
img_quarnt=[]
quant_levels = [2, 4, 8, 16, 32, 64, 128]
for q_level in quant_levels:  
    img_quarnt.append(quantize(img, q_level))


fig,ax = plt.subplots(2,4,figsize=(15,15))
ax[0,0].set_title("orig")
ax[0,1].set_title("quant level = 128")
ax[0,2].set_title("quant level = 64")
ax[0,3].set_title("quant level = 32")
ax[1,0].set_title("quant level = 16")
ax[1,1].set_title("quant level = 8")
ax[1,2].set_title("quant level = 4")
ax[1,3].set_title("quant level = 2")

ax[0,0].imshow(img, cmap = 'gray', vmin = 0, vmax = 255,aspect='equal')
ax[0,1].imshow(img_quarnt[6], cmap = 'gray', vmin = 0, vmax = 255,aspect='equal')
ax[0,2].imshow(img_quarnt[5], cmap = 'gray', vmin = 0, vmax = 255,aspect='equal')
ax[0,3].imshow(img_quarnt[4], cmap = 'gray', vmin = 0, vmax = 255,aspect='equal')
ax[1,0].imshow(img_quarnt[3], cmap = 'gray', vmin = 0, vmax = 255,aspect='equal')
ax[1,1].imshow(img_quarnt[2], cmap = 'gray', vmin = 0, vmax = 255,aspect='equal')
ax[1,2].imshow(img_quarnt[1], cmap = 'gray', vmin = 0, vmax = 255,aspect='equal')
ax[1,3].imshow(img_quarnt[0], cmap = 'gray', vmin = 0, vmax = 255,aspect='equal')
plt.tight_layout()
plt.show()


#%%



fig,ax = plt.subplots(2,4,figsize=(15,10))
ax[0,0].set_title("orig")
ax[0,1].set_title("quant level = 128")
ax[0,2].set_title("quant level = 64")
ax[0,3].set_title("quant level = 32")
ax[1,0].set_title("quant level = 16")
ax[1,1].set_title("quant level = 8")
ax[1,2].set_title("quant level = 4")
ax[1,3].set_title("quant level = 2")

ax[0,0].hist(img.flatten(), 256, [0, 256])
ax[0,1].hist(img_quarnt[6].flatten(), 256, [0, 256])
ax[0,2].hist(img_quarnt[5].flatten(), 256, [0, 256])
ax[0,3].hist(img_quarnt[4].flatten(), 256, [0, 256])
ax[1,0].hist(img_quarnt[3].flatten(), 256, [0, 256])
ax[1,1].hist(img_quarnt[2].flatten(), 256, [0, 256])
ax[1,2].hist(img_quarnt[1].flatten(), 256, [0, 256])
ax[1,3].hist(img_quarnt[0].flatten(), 256, [0, 256])
fig.tight_layout()
plt.show()

Equation of decimation and interpolation to the 3 different methods:

sinc:

2.2.2-3 in matlab